
clear all;
clc;
disp('forecast one-quarter ahead');
load optimization_g1000k.mat

INSAMPLE_YR = 2008;
INSAMPLE_QR = 4;
% CRITERIA = 10;

inSampleLength = (INSAMPLE_YR - 1998)* 4 + INSAMPLE_QR;
outSampleLength = NUMBER_OF_YEAR * Quarter_LIST_LENGTH - inSampleLength;


% xx =zeros(RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
% avgRollingTM =zeros(length(SVD),RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
% bin=zeros(length(SVD),RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
% InsampleSVD = zeros(inSampleLength,1);
% InsampleMobilityNorm = zeros(inSampleLength,1);
% PHat=zeros(10,inSampleLength,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
% N = zeros(inSampleLength+outSampleLength,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
% naive1 = zeros(inSampleLength, RATE_LIST_LENGTH-1, RATE_LIST_LENGTH );
% naive2 = zeros(inSampleLength, RATE_LIST_LENGTH-1, RATE_LIST_LENGTH );

filename='optimization_08q4_jun13.xls';

% To calculate credit cycle index Z, real and forecast----------------------
% To minimize the distance to get estimate of w----------------------------

% Website tutoring upload data from xls to matlab:
% http://blinkdagger.com/matlab/matlab-using-xlsread-to-import-excel-data/

% [fcstSVD] = xlsread('fittedReg_SampleEndTo_08Q4_09Q3_SVD.xlsx','B50:B53');
[fcstInvPD] = xlsread('fittedReg_SampleEndTo_08Q4_09Q3_PD.xlsx','K51:K54');

% PHat_BFS = zeros (inSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);


P = stdProbQ;
Z_fcst = zeros (outSampleLength, 1);
Wfcst = zeros (outSampleLength, 1);
Pfcst = zeros(outSampleLength,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
efcst = zeros(outSampleLength,CRITERIA,3);
g = zeros(outSampleLength,4);
s = zeros(outSampleLength,4);

for i = 1: outSampleLength ;
    % To calculate credit cycle index Z, real and forecast
    %----------------------------------------------------------% D3
        PD = defaultFreq(1:inSampleLength+i-1);
        inversePD = norminv(PD);
        Z = - zscore(inversePD);
        Z_fcst(i) =  - ( fcstInvPD(i) - mean(inversePD) ) / std(inversePD);
    %----------------------------------------------------------% SVD
%     SVD_temp = SVD(1:inSampleLength+i-1);
%     Z = - zscore(SVD);
%     Z_fcst(i) =  - ( fcstSVD(i) - mean(SVD_temp) ) / std(SVD_temp);
%     
%     t = inSampleLength+i-1;
%     X(:,:) = bin(inSampleLength+i-1,:,:);
    %------------------------------------------------------------
    
    % To minimize the distance to get estimate of w
    obj = @(wfcst)difD3_fixedW(wfcst,t,P,X,Z,R,C); % D3
%     obj = @(wfcst)difSVD_fixedW(wfcst,t,P,X,Z,R,C); % SVD
    [foundW,fval,exitflag,output] = fminbnd(obj,0,1);
    Wfcst(i,1) = foundW;
    % subplot(5,2,8); fplot(obj8,[0 1],'*');
    % title('Mininization Object with D3 distance',...
    %   'FontWeight','bold');
    
    %To calculate out-of-sample fcst matrix
    Z(inSampleLength+i) = Z_fcst(i);
    Ptemp= calculatePHat_fixedW(Wfcst(i,1),t,X,Z,R,C); % PHat = f(X,What,Zhat)
    Pfcst(i,:,:) = Ptemp (t,:,:);
    clear Ptemp
    
    
    % error(Mean Absolute Error)
    
    
    x1(:,:) = Pfcst(i,:,:);
    x2(:,:) = avgRollingTM (inSampleLength+i-1,:,:); % take average as naive1 benchmark
    x3(:,:) = stdProbQ (inSampleLength+i-1,:,:);% take previous as naive2 benchmark
    y(:,:) = stdProbQ(inSampleLength+i,:,:);
    for type = 1: CRITERIA
        efcst(i,type,:) = calculateMAE(type,x1,x2,x3,y,R,C);
    end
    
    % visual check est. TM
    %---------------------------------------------------------------D3
    for ii = RATE_LIST_LENGTH-1
        for jj = RATE_LIST_LENGTH
    s(i,1) = s(i,1) + (ii - jj)* sign(x1(ii,jj) - y(ii,jj))*(x1(ii,jj) - y(ii,jj))^2;
    s(i,2) = s(i,1) + (ii - jj)* sign(x2(ii,jj) - y(ii,jj))*(x1(ii,jj) - y(ii,jj))^2;
    s(i,3) = s(i,1) + (ii - jj)* sign(x3(ii,jj) - y(ii,jj))*(x1(ii,jj) - y(ii,jj))^2;
        end
    end
%-------------------------------------------------------------------SVD
%     s (i,1) = mean(svd(x1)); s(i,2) = mean(svd(x2)); s(i,3) = mean(svd(x3)); s(i,4) = mean(svd(y));
%---------------------------------------------------------------------------     
     
    g (i,1)= mean(x1(:,C));  g(i,2) = mean(x2(:,C)); g(i,3) = mean(x3(:,C)); g(i,4)= mean(y(:,C));% average default
    gQ (:,1) = x1(:,C); gQ (:,2) = x2(:,C); gQ (:,3) = x3(:,C); gQ (:,4) = y(:,C); %each rate default in every quarter
    
    subplot (2,2,i);plot (gQ);
    legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');
    title('Default Probability Forecast: 2009Q1 - 2009Q4');
    
end

figure
plot (g);
legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');
set(gca,'XTick',1:1:4);
set(gca, 'XTickLabel', {'09Q1','09Q2','09Q3','09Q4'});
title('Default Probability Forecast: 2009Q1 - 2009Q4');


figure
plot (s);
legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');
set(gca,'XTick',1:1:4);
set(gca, 'XTickLabel', {'09Q1','09Q2','09Q3','09Q4'});
title('SVD Compariaon: Forecast and Real: 2009Q1 - 2009Q4');

% % visual check est. TM-----------------------------------------------------
%
% % first, we look at criteria D3
% g1 (:,1)= Pfcst(i,1,1);  g1(:,2) = naive1(:,1,1); g1(:,3) = naive2(:,1,1); g1(:,4)= stdProbQ (1:length(PHat),1,1);% rank 1.5 - 1.5
% g2 (:,1)= PHat(8,:,3,3); g2(:,2) = naive1(:,3,3); g2(:,3) = naive2(:,3,3); g2(:,4)= stdProbQ (1:length(PHat),3,3);% rank 2.5 - 2.5
% g3 (:,1)= PHat(8,:,6,6); g3(:,2) = naive1(:,6,6); g3(:,3) = naive2(:,6,6); g3(:,4)= stdProbQ (1:length(PHat),6,6);% rank 3 - 3
% g4 (:,1)= PHat(8,:,5,8); g4(:,2) = naive1(:,5,8); g4(:,3) = naive2(:,5,8); g4(:,4)= stdProbQ (1:length(PHat),5,8);% rank 3.5 -D
% g5 (:,1)= PHat(8,:,6,8); g5(:,2) = naive1(:,6,8); g5(:,3) = naive2(:,6,8); g5(:,4)= stdProbQ (1:length(PHat),6,8);% rank 4 -D
% g6 (:,1)= PHat(8,:,7,8); g6(:,2) = naive1(:,7,8); g6(:,3) = naive2(:,7,8); g6(:,4)= stdProbQ (1:length(PHat),7,8);% rank 4.5 -D
%
% % then, we look at criteria SVD
%
% g6 (:,1)= PHat(5,:,7,8); g6(:,2) = naive1(:,7,8); g6(:,3) = naive2(:,7,8); g6(:,4)= stdProbQ (1:length(PHat),7,8);% rank 4.5 -D
%
% figure
% plot (g6); %([g1,gg1,ggg1]);%,'-*'); plot (gg1,'-rs');
% set(gca,'XTick',1:4:length(g1));
% labels = quaterlabels(1998, length(g1));
% labels = labels(1:4:length(g1));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% legend('model','benchmark1-avg','benchmark2-previous','empirical','Location','NorthWest');




% % W = zeros (outSampleLength, 3); % first column is for metric BFS, the second colum is for metric D2
% PHat_BFS = zeros (inSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
% % PHat_D2 = zeros (outSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
% % PHat_SVD = zeros (outSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
% R = RATE_LIST_LENGTH-1;
% C = RATE_LIST_LENGTH;
% yrEst= INSAMPLE_YR - 1998 + 1 + floor(INSAMPLE_QR / 4);
%
% % for i = 1: inSampleLength %outSampleLength
%
% P = stdProbQ;
% X(:,:) = bin(inSampleLength,:,:);
% wantedZ = Z;
% t = inSampleLength;
%
% % BFS distance
% obj1=@(w)difBFS_fixedW(w,t,N,P,X,wantedZ,R,C);
% [foundW,fval,exitflag,output] = fminbnd(obj1,0,0.5);
% figure;
% fplot(obj1,[-0.5 0.5],'*')
% title('Mininization Object with BFS distance',...
%   'FontWeight','bold')
% PHat_BFS = calculatePHat_fixedW(foundW,t,X,wantedZ,R,C);
%
%
% % W(i,1) = foundW;
% % sheetlist=strcat('',year_list_text(yrEst,:),'.',quarter_list_text(qrEst,:),'');
% % xlswrite(filename, PHat, sheetlist);
% % figure;
% % subplot(outSampleLength,1,i); fplot(obj1,[-0.5 0.5],'*')
% % title('Mininization Object with BFS distance',...
% %   'FontWeight','bold')
%
%
% % D2 distance
% obj2=@(w)difD2(w,P,X,wantedZ,R,C);
% [foundW,fval,exitflag,output] = fminbnd(obj2,0,0.5);
% W(i,2) = foundW;
% PHat_D2(i,:,:) = calculatePHat(foundW,X,wantedZ,R,C);
%
%
% % SVD
% obj3=@(w)difSVD_fixedW(w,t,P,X,wantedZ,R,C);
% [foundW,fval,exitflag,output] = fminbnd(obj3,0,0.5);
% % W(i,3) = foundW;
% PHat_SVD = calculatePHat_fixedW(foundW,t,X,wantedZ,R,C);
% % xx(:,:) = PHat_SVD(i,:,:);
% % sheetlist=strcat('',year_list_text(yrEst,:),'.',quarter_list_text(qrEst,:),'');
% % xlswrite(filename, xx, sheetlist);
% % clear xx;
% % end




% MAE ---------------------------------------------------------------------
% mae_perf = zeros(5,outSampleLength);
% mae_svd_perf = zeros(5,outSampleLength);
%
% for i = 1: outSampleLength
%     P(:,:) = stdProbQ(inSampleLength+i,:,:);
%     previous (:,:) = stdProbQ(inSampleLength+i-1,:,:);
%     average (:,:) = avgRollingTM(inSampleLength+i-1,:,:);
%     model_BFS (:,:) = PHat_BFS(i,:,:);
%     model_D2 (:,:) = PHat_D2(i,:,:);
%     model_SVD (:,:) = PHat_SVD(i,:,:);
%
%
%     mae_perf(1,i) = mean2 (abs (P - previous));
%     mae_perf(2,i) = mae (P - average);
%     mae_perf(3,i) = mae (P - model_BFS);
%     mae_perf(4,i) = mae (P - model_D2);
%     mae_perf(5,i) = mae (P - model_SVD);
%
%     mae_svd_perf(1,i) =  abs (mean(svd(P))- mean(svd(previous)));
%     mae_svd_perf(2,i) =  abs (mean(svd(P))- mean(svd(average)));
%     mae_svd_perf(3,i) =  abs (mean(svd(P))- mean(svd(model_BFS)));
%     mae_svd_perf(4,i) =  abs (mean(svd(P))- mean(svd(model_D2)));
%     mae_svd_perf(5,i) =  abs (mean(svd(P))- mean(svd(model_SVD)));
%
%
% end

% naive1=stdProbAvgQ;
% naive2=stdProbQ(inSampleLength+i,:,:);
% % n(:,:)=N(yrp,freqjp,2:10,13);
% % nn=repmat(n,1,10);
% % nt(1,1,:,:)=nn;
% % Preal=P(:,:,2:10,2:11);
% e_mod=(nt.*(Preal-PZ).^2)./(PZ.*(I-PZ));
% e_avg=(nt.*(Preal-naive1).^2)./(naive1.*(I-naive1));
% e_pre=(nt.*(Preal-naive2).^2)./(naive2.*(I-naive2));
% perf_mod = mae(e_mod);
% perf_naive1 = mae(e_avg);
% perf_naive2 = mae(e_pre);
% %
%
% Parameters={'foundW','Zhat', 'perf_mod', 'perf_naive1','perf_naive2'; foundW Zhat perf_mod perf_naive1 perf_naive2}';
% xlswrite(filename, Parameters,'Parameters');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Naive1 MAE of L1,L2 NSD,D1,D2   %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% L1_n1=Preal-naive1;
% L2_n1=(Preal-naive1).^2;
% NSD_n1=((Preal-naive1).^2)./(naive1);
% % D1_n1=
% % D2_n1=
%
% perf_naive1e_L1=mae(L1_n1);
% perf_naive1e_L2=sqrt(mae(L2_n1).*90)/90;
% perf_naive1e_NSD=mae(NSD_n1);
% % perf_naive1e_D1=
% % perf_naive1e_D2=
%
% disp('fini');




% for i = 1 : outSampleLength
%
% % SVDhat = SVD (1:inSampleLength + i -1, :);
% % SVDhat(inSampleLength + i, :) = fcstSVD (i, :);
% % ZhatVector = - zscore (SVDhat);
% % Zhat(i, :) = ZhatVector(size(ZhatVector,1),:);
% % save in 'optimization_07q4_jun2_3.xls';
%
% Zhat(i,:) = - ( fcstInvPD(i,:)- mean(inversePD(1:inSampleLength+i-1,:)) )/std(inversePD(1:inSampleLength+i-1,:));
% % save in 'optimization_07q4_jun2_5.xls';
%
% end

